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ABSTRACT 

Analytic estimates of the viscous time-scale due to cloud-cloud collisions have 
been as high as thousands of Gyr. Consequently, cloud collisions are widely ignored as 
a source of viscosity in galactic disks. However, capturing the hydrodynamics of discs 
in simple analytic models is a challenge, both because of the wide dynamic range and 
importance of 2D and 3D effects. To test the validity of analytic models we present 
estimates for the viscous time-scale that are measured from three dimensional SPH 
simulations of disc formation and evolution. We have deliberately removed uncertain- 
ties associated with star-formation and feedback thereby enabling us to place lower 
bounds on the time-scale for this process. We also contrast collapse simulations with 
results from simulations of initially stable discs and examine the impact of numerical 
parameters and assumptions on our work, to constrain possible systematics in our 
estimates. We find that cloud-collision viscous time-scales are in the range of 0.6-16 
Gyr, considerably shorter than previously estimated. This large discrepency can be 
understood in terms of how the efficiency of collisions is included in the analytical 
estimates. We find that the viscous time-scale only depends weakly on the number 
of clouds formed, and so while the viscous time-scale will increase with increasing 
resolution, this effect is too weak to alter our conclusions. 
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1 INTRODUCTION 

Arguably, the most successful model for the formation of disc 
galaxies is the ACDM model, in which galaxies are formed 
from the dissipational collapse of baryonic gas within a dark 



by feeding this turbulence (Vollmer & Beckert 20031. Ad- 
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the physical viscosity of the baryonic gas is not anticipated 
to have a strong influence on gas evolution except in magne- 



tized or hot environments such as a galaxy cluster ( Sijacki & 



Springel |2006[ ), effective kinematic viscosities could in prin- 
ciple impact disc evolution. Simulations by |Lin fc Pringle| 
( 1987 1 with a viscous time-scale close to the star-formation 



time-scale showed that viscous evolution with infall can re- 
produce the ubiquitous exponential density profile from a 
range of initial conditions. In this work the viscosity was as- 
sumed to be caused by large-scale turbulent motions dissi- 
pating kinetic energy and transporting angular momentum. 
Feedback from supernovae can be a source of viscosity 
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ditionally, the self-gravity of the gaseous disc can provide 

This can 



an effective viscosity ( Vollmer & Beckert 2002 1 
take the form of large-scale instabilities ( Rafikov 2009 , Gam 
|2001 1, or of interactions between giant molecular clouds 



} Vollmer &c Beckert||2002 1. These cloud interactions poten- 
tially generate viscosity through two different mechanisms 
Firstly, gravitational scattering can increase the velocity dis- 
persion of the cloud population, converting orbital energy 
into large-scale turbulence ( Gammie et al.|[l991 Fukunaga 



fc Tosa|1989| |Agertz et al.|2009[ ). Secondly, during inelastic 



collisions between clouds, shocks convert orbital energy into 
turbulence and heat within the colliding clouds (e.g. |Git- 



|tins et al.|2 003 Kitsion as~fe Whitw orth 2007, Anathpindika 
2009). Radiative processes contribute to the dissipation of 
kinetic energy during these collisions, and are also impor- 
tant for dissipating turbulent energy that has cascaded into 
thermal energy. These processes are significant even in the 
absence of star-formation: the observations compiled by |Dib| 
|et al.| ( |2006[ ) show that the velocity dispersion of HI gas does 
not strongly depend on the star-formation rate below a cer- 
tain threshold, and the AMR simulations of |Agertz et al.| 
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( 2009 1 suggest that a 'baseline' turbulence is caused by inter- 



actions between clouds, and that this is only supplemented 
by supernova feedback at high star-formation rates. 



It has been argued (Vollmer & Beckert 2002 Bell 



2002 hereafter B02) that cloud-collisions are not an efficient 



source of viscosity. In particular, in B02 the time-scales for 
viscosity due to cloud collisions are estimated to be on the 
order of t v ~ 1000 Gyr in most local spiral galaxies, although 
the time-scales might be considerably lower in earlier gas- 
rich galaxies or in galaxies where the velocity distribution 
of GMCs has been stirred up by some mechanism (such as 



galaxy interaction e.g. Hernquist & Mihos 1995 I. Vollmer & 



Beckert ( 2002 1 argues that because molecular clouds evap- 

10 7 yr, and this is less than the time 

10 8 yr), cloud collisions are very rare. 



orate at an age of <- 
between collisions {>• 
However, cloud formation times, assuming that collapse and 
formation of H2 are the dominant factors in forming a cloud, 



appear to equally short (Glover & Mac Low 20071. This 



leads to a scenario in which the number density of clouds 
is roughly constant, although the short life-time may affect 
the velocity dispersions of molecular clouds as they have 
less time to build up a large deviation from circular velocity 
through scattering events with other clouds. In this steady 
state, the effective collision time-scale should remain similar. 

It has also been argued that physical collisions between 
clouds have a smaller effect than gravitational scattering 
I Jog & Ostrikcr 1988), particularly as if magnetic fields are 
taken into account, cloud interaction cross-sections may be 



underestimated (Ozernoy et al. 19981. On the other hand, 



Das & Jog ( 1996 1 modelled a system of cloud particles, find- 
ing that cloud collisions rather than local gravitational inter- 
actions (scattering events) dominate the mass distribution 
and velocity dispersion of molecular clouds, suggesting that 
cloud collisions may indeed be important. However, as far 
as we are aware, the effective viscosity of direct cloud-cloud 
collisions has not yet been examined in global three dimen- 
sional numerical hydrodynamic models. 

Most simulations of cloud formation and the associ- 
ated disc dynamics have been performed in two dimensions 



and/or on a small scale using shearing-box studies {e.g. Kim 
fc Ostriker|2007 1. However, increased computing power and 



the availability of locally adaptive algorithms have recently 
enabled galaxy-scale simulations with sufficiently high reso- 
lution to resolve cloud-formation in discs. Numerical experi- 



ments have been performed using both AMR ( Tasker & Tan 
2009{|Tasker|2011| |Agertz et al.|2009[ ) and SPH ( jRobertson 
fc Kravtsov|2008 1 with resolutions as fine as 6 pc. The non- 



trivial cooling processes and chemistry make these simula- 
tions a significant technical challenge. Agertz et al. ( 2009 1 



and Tasker & Tan ( 2009 1 ran suites of high resolution AMR 



simulations of Milky- Way and M33-like disc galaxies, and 
reported on the properties of the clouds generated by their 
models, including cloud-cloud velocity dispersion. However, 
neither study has provided an estimate of the viscous time- 
scale due to cloud-cloud collisions. Furthermore, the discs of 



Tasker & Tan ( 2009 1 are much more stable than the Milky 



Way, with a density distribution chosen to give a constant 
value of the Toomre Q parameter (Toomre 19641, and a 



static dark matter and stellar component, which may in- 
hibit some of the instabilities important to cloud formation. 

In this paper we revisit the calculations of B02 with 
full three dimensional SPH models. This is not entirely triv- 



ial since there is no universally agreed upon cloud finding 
process. However, the use of a particle method enables the 
Friends-of- Friends ( Davis et al.|1985 1 group finding method- 
ology and we adapt that to our simulations. Hence given 
our cloud population our primary goal is to see whether the 
analytical calculations are supported, and if not what the 
implications are. It is important to note that the results of 
such simulations could highlight non-physical evolution in 
numerical schemes with artificial viscosities, of which SPH 
is a notable example ( Valdarnini|2011 1 . We also investigate 
the issue of numerical artefacts in our calculated results. 
This is a key issue since structure formed within simula- 
tions starting from smooth initial conditions is inevitably 
the result of amplification of noise in the initial conditions. 

While a full calculation in the cosmological context {e.g. 
Katz 1992; Thacker & Couchman 200T1 IGovernato et al. 



|2007[ |Stinson et al7 | |2010 | |Scannapieco et al.||2009) [Brook 



|et al.|2008 1 is beyond the scope of this paper, primarily due 
to resolution limitations, we instead consider two classes of 
isolated models. We examine an equilibrium system with 
similar parameters to the Milky Way consisting of a gas disc, 
a stellar disc and bulge, and a dark matter halo. Here the 
gas disc is stabilized by the other components which dom- 
inate the system's mass. We also consider the dissipational 
collapse problem that has been used extensively elsewhere 
{e.g. |Gott & Thuan||1976| |Carlberg||1984| |Katz & Gunn 



1991 Brook et al. 2004 Kaufmann et al. 20061. By con 



trast with the Milky Way model, this collapse produces a 
very unstable disc, and so we investigate both high-stability 
gas-poor systems and low-stability gas-rich systems. These 
models include hydrodynamics, gravitational interactions, 
and cooling with a dynamic temperature floor. By removing 
the numerous unknowns associated with star-formation and 



feedback (as discussed in numerous places e.g. Thacker & 



Couchman|2000||Ceverino fc K lypin 2009; Christcnsc n et al 



2010 1 we hope to isolate the impact of cloud-cloud interac- 



tions and place lower bounds on the viscous time-scale. 

The layout of the paper is as follows: in section [2] we 
outline the details of our simulation code. We also discuss 
the initial conditions, our cloud identifying approach and 
also the underlying theory of the effective viscosity. Results 
are presented in section [3] followed by a brief conclusion. 



2 SIMULATION 



2.1 Simulation Code 



We model the dark matter, stars and gas using a specially 
adapted version of the OpenMP n-body AP 3 M ( Couchman 
T99"Tj) SPH (|Monaghan||T992"| > code HYDRA ( jThacker fc 
Couchman 20061. Our modifications are: 



(i) The cooling curve has been extended to down to 10 K 



using the cooling function (A) of Wada & Norman ( 2001 1, al- 
though we set our fiducial temperature floor to 300K to make 



our results more comparable with Tasker & Tan (2009), ex- 



cept in cases where we investigate the effect of a lower floor. 



The earlier cooling curve of Sutherland & Dopita ( 1993 1 is 
retained for T > 10 4 K. The combined cooling curve is plot- 
ted in Fig. [1] 

(ii) We implemented a dynamic temperature floor based 
on Robertson & Kravtsov (20081, described in section 2.1.1 
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2.2 Initial Conditions 
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Figure 1. The cooling curve used in ou r models. Values below 
10 4 K are from|Wada fc Norman" ||200l| , while those above 10 4 
K are from [Sutherland fc Dopita| 1993^ 



(iii) The parallelisation algorithm has been altered so that 
during the particle-particle gravity and SPH calculation re- 
gions containing a large number of particles are split over 
all processors, instead of each processor receiving a single 
region. This greatly improves load-balance in simulations 
containing many dense clumps of particles. 



2.1.1 Dynamic Temperature Floor 



We use a method similar to Robertson & Kravtsov ( 2008 1 



to ensure that the Jeans mass is resolved in our simulations. 



This is to satisfy the Truelove et al. ( 1997 1 criterion and 



avoid artificial fragmentation - crucial in simulations of cloud 
formation. The method is in the form of a dynamic pressure 
floor. The Jeans mass( Jeans|1 902) is defined as 



6G 3/ 2p i/2 ■ 



(1) 



Bate & Burkert ( 1997 1 noted each particle should sat- 
isfy 2N neigh m gas < m Jeans (where N neigh is the number of 
SPH neighbours for the particle and m gas is the gas parti- 
cle mass) to avoid artifical fragmentation. Defining the local 
ratio of the Jeans mass to the SPH kernel mass as /ij ea ns, 
then this requirement is equivalent to the statement that 



h 3 



6G 3 / 2 N n 



9 X/2 



< Nj, 



(2) 



* neighTTlgas P 

with A^jeans being the required factor by which the Jeans 
mass must be resolved, which in the Bate and Burkert case 
is set to 2. In an ideal gas c s oc y/u, so we can fulfill this 
criterion by applying 

2/3 



(3) 



whenever /ij ca ns < Nj ea . ns - 

We found spurious (sub-resolution) string-like struc- 
tures forming within clouds for low values of -/Vj cans , and 
found that iV,j cans = 50 removed these structures and re- 
sulted in a more homogeneous interior for clouds. 



2.2.1 Milky Way Model 

We produce our Milky Way model using the GALACTICS 



package (Kuijken & Dubinski 1995 Widrow & Dubinski 
|2005| |Widrow et al.||2008[ > with the parameters in Table 2 
of |Widrow et al.| ( |2008| ) .Through an iterative process, this 
package produces an equilibrium system consisting of an ex- 
ponential stellar disc, a stellar bulge and a dark matter halo. 
The disc is exponential radially, follows sech 2 vertically, and 
has a radial dispersion profile of 



o-r(R) = o- 2 R0 exp(-i?/i? CT ), 



(4) 



where we set R a — Rd for simplicity. We generate the gas 
disc by copying the disc star particle positions and flipping 
the coordinates across the x — y plane to prevent particles 
having coincident positions. Bulge particles are not copied. 
The masses of the gas and star particles are scaled to give 
the appropriate mass ratio. The gas disc is given a disper- 
sionless velocity profile output by GalactlCs and is initially 
isothermal at 10 4 K. The disc scale length is 2.81 kpc, trun- 
cated at 30 kpc with a truncation scale-length of 0.1 kpc. 
The scale height is initially 0.36 kpc, and the total disc mass 
is 5.2 x 1O 1O M . The halo density profile is 



P = 



'01 



{r/a h )"<(l + r/a h y-i 2 



erfc 



rh 



(5) 



\V2Sr h , 

where is the halo scale parameter, r^ is the cutoff radius, 
Srn is the scale length for the truncation, 7 is the 'cuspiness' 
parameter (equal to unity for an NFW profile), and is a 
velocity parameter that sets the mass of the halo. 

Halo parameters are given in Table [I] We also ran a 
test with a static analytic halo potential, to explore if the 
discretization of the halo has any effect on cloud formation. 

The bulge density profile is 



Pb 



(6) 



where p = 1 — 0.6097/n + 0.05/n 2 gives a Sersic profile with 
n the Sersic index. R e is the radial scale parameter, and in 
GALACTICS pb is parametrized by the velocity parameter 
a b = {47rn6 n(p - 2) r[7i(2 ~ p)]R 2 e p b } 1/2 . We set these param- 
eters to n = 1.31, at — 272 km s , and R e = 0.64 kpc. 

We have named our fiducial run LowSoftMW. To test 
the effects of change the resolution, softening length, tem- 
perature floor, gas mass fraction, and artificial viscosity we 
investigate a total of ten different runs, summarized in Ta- 
ble[2] Both MidSoftMW and HighSoftMW have higher grav- 
itational softening lengths; MedGasMW and HighGasMW 
have higher gas mass fractions; LowFloorMW has a lower 
temperature floor; LowViscMW has lower artificial viscosity 
parameters (a, /3); and LowResMW has a lower resolution. 
In addition, as a convergence check we ran a higher resolu- 
tion simulation (HighResFlatMW) with a total of 3.5 x 10 6 
particles and a softening length of 45 pc, although we do 
not consider this our fiducial run as the simulation time 
did not reach a full Gyr. We found when running a simula- 
tion of this high resolution with identical initial conditions 
to LowSoftMW that the disc was dominated by a strong 
ring-shaped shock propagating outwards. This is caused by 
a combination of the rapid vertical collapse of the disc as it 
initially cools, and the rotation curve not being quite precise 
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Name 


a h (kpc) 


Th (kpc) 


Sr h (kpc) 


7 


a h (km s 1 ) 


Mhalo 


Collapse haloes 


25.75 


300 


50 


1.0 


351 


1.1 X 10 12 


MW haloes 


13.6 


275 


25 


0.81 


330 


7.3 x 10 11 



Table 1. Halo parameters. As in Eq. [3] is the halo scale parameter, is the truncation radius, Sr^ is the scale length for this 
truncation, 7 is the cuspiness parameter, and is a velocity parameter that sets the halo mass, Mh a i . 



enough as GALACTICS is intended for collisionless mechan- 
ics and does not take into account gravitational softening or 
the pressure gradient of the gaseous disc. At the lower reso- 
lutions this shock is not well captured, and the disc quickly 
returns to equilibrium, so this is only a problem at our high- 
est resolution. To prevent the shock becoming a problem it 
is necessary to start the simulation from an initially flat- 
tened state akin to the later evolution of the cooled disks. 
We therefore flattened the gas disc to a similar scale height 
as the cooled disks, which is a factor of 10 smaller. Circular 
velocities (w c irc) were then set up using radial accelerations 
((Jrad) generated from a single iteration of the HYDRA code, 
and explicitly setting a ra d = v^ ilc /R for each gas particle, 
where R is the radial coordinate of the particle. We also 
performed a simulation (FlatMW) with these initial condi- 
tions but at our fiducial (moderate) resolution, for a fair 
comparison of the effects of resolution. 



2.2.2 Monolithic Collapse Model 

This model consists of a spherically symmetric distribution 
of gas within an equilibrium NFW dark matter halo. We 
generate the halo using GALACTICS according to the pa- 

''Ma. 



rameters in Table [T] giving a halo with M = 1.1 X 10 1 

For the gas we use the 'high-entropy' profile of Kauf- 
mann et al. ( 2009|, which was produced from equation 1 
of Kazantzidis et al. (20041, setting c = 1, a = 1, ft = 3, 



and 7 = 0. Kaufmann noted that a gas density profile that 
is shallower than the NFW profile (as expected in models 
with pre-heating feedback e.g. Mo & Mao 2002) produces 



an angular momentum distribution in the final object that 
better fits observations. In this model, the gas collapses into 
clumps which combine to form an unstable disc. 

As in Kaufmann et al. (20071, the initial temperature 



profile is calculated to provide hydrostatic equilibrium ac- 
cording to 



T(r) = 



/' 



1 



J . .GMtot(r) , 
P^ r ) ~2 dr > 



(7) 



where [i is the mean molecular weight of the gas (taken as its 
primordial value, fi ~ 0.59m_f/), ks is the Boltzmann con- 
stant, pG is the initial gas density, and M to t(r) is the total 
mass (gas and dark matter) within a sphere of radius r. We 
give the gas a flat velocity profile. The positions of the gas 
particles in our initial conditions are simply the generated 
positions of the dark matter particles flipped as in section 

Em 

To set up a rotating halo, GALACTICS swaps a frac- 
tion of the dark matter particles' velocities over the radial 
axis to increase the number of particles rotating in the same 
direction. We assume the gas and the dark matter have the 
same specific rotational momentum, i.e. 



\Lg\ 



(8) 



Mdm 

so that the spin parameter ( Binney & Trcmainc 2008 ) of the 
gas is equal to the spin parameter of the halo, 



\Lq\ \/|-Edm| _ I I/dm I \/\Emj[\ 



Mg GM 3/2 



M DM GM^ 



Ar 



(9) 



We used a spin parameter of \g = 0.038, close to the median 
value observed in simulations ( |Bullock et al. 2001 Barnes 
& Efstathiou 19871. After each gas/DM halo is produced, 



it is evolved for 0.5 Gyr with cooling switched off to ensure 
the ICs are stable. Our first model (HighSoftC) is performed 
with the softening equal to Kaufmann et al. |2009|'s and the 



temperature floor equal to Kaufmann et al.| (|2009j)'s cool- 
ing floor. We also investigated models with lower softening 
lengths and temperature floors to see if smaller clouds were 
resolved. A low resolution run was performed as a conver- 
gence check, and finally we performed a model with a low 
gas fraction to see the effect of increasing the disc's stability. 
These models are summarised in Table [3] 



2.3 Cloud Identification & Analysis 

2.3.1 Identification algorithm 

We identified clouds using a two-step process. A density 
threshold was applied and then particles above this thresh- 
old were linked into groups using the friend-of-friends algo- 
rithm. Using a density threshold partially avoids the noto- 
rious 'string of pearls' effect that that may lead to spurious 
filamentary structures or the merging of many smaller struc- 
tures into one larger one. We set the density threshold at a 
density of n = 7 Mqpc -3 , and used a linking length of 50 
pc. 

The high density threshold ensures that only dense 
cloud-like objects are selected, while the linking length is 
close to the size of the softening ensuring small fluctuations 
below this threshold can be skipped over. We set the mini- 
mum cloud size to 30 particles, giving minimum cloud mass 
of 1.6 x 10 5 -8.5x 1O 6 M , depending on resolution. We found 
these parameters largely selected dense, cool clumps while 
excluding other objects such as filaments from cloud encoun- 
ters. Because our lower limit for cloud size is a number of 
particles and not a mass, we include clouds of increasingly 
small mass as we increase the resolution of our simulation. 
This may create a resolution dependence until individual 
molecular clouds are resolved - a level that we do not achieve 
in this work, although as noted in section |3.1.2| the major 
axisymmetric instabilities appear to be resolved. 

Clouds are tracked from output to output by examining 
the particles resident in each cloud. If the cloud A at time ti 
contains at least half of the particles contained by cloud B 
at the time of the following output U+i, then A is a parent 
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I rt !oc) 

t sojt / 


Tr (K) 

* joor V / 


71* 




"dm 


772g j TYh* 


t a (Gvr) 


"-disc V i / 


(X 


8 


LowSoftMW 


60 


300 


5 x 10 5 


4 x 10 5 


5 x 10 5 


0.1 


1.116 


0.36 


1, 


2 


Medbott M W 


100 


300 


5 x 10 J 


4 x 10 5 


5 x 10 5 


0.1 


1.146 


0.36 


1. 


2 


HighSoftMW 


500 


300 


5 x 10 5 


4 x 10 5 


5 x 10 5 


0.1 


1.542 


0.36 


1. 


2 


LowResMW 


60 


300 


1 x 10 5 


8 x 10 4 


1 x 10 5 


0.1 


1.959 


0.36 


1, 


2 


LowFloorMW 


60 


10 


5 x 10 5 


4 x 10 5 


5 x 10 5 


0.1 


1.004 


0.36 


1, 


2 


LowViscMW 


60 


300 


5 x 10 5 


4 x 10 5 


5 x 10 5 


0.1 


1.002 


0.36 


0.5 


. 1 


MedGasMW 


60 


300 


5 x 10 5 


4 x 10 5 


5 x 10 5 


0.2 


0.485 


0.36 


1, 


2 


HighGasMW 


60 


300 


5 x 10 5 


4 x 10 5 


5 x 10 5 


0.5 


0.434 


0.36 


1, 


2 


FlatMW 


60 


300 


5 x 10 5 


4 x 10 5 


5 x 10 5 


0.1 


0.790 


0.036 


1, 


2 


HighResFlatMW 


45 


300 


1.25 x 10 6 


1 x 10 6 


1.25 x 10 6 


0.1 


0.318 


0.036 


1, 


2 



Table 2. Summary of Milky Way runs. l 30 ft is the minimum softening length, Tfi oor is the temperature floor, n„, n g , and hdm are the 
numbers of star, gas and dark matter particles, m g /m* is the gas/star mass ratio for the disc, t enl i is the total simulation time, /i d i sc is 
the scale height of the disc, and a and 8 are the artificial viscosity parameters. 



Name 


hoft (pc) 




DOT (-^) 




n g 


"DM 


m g /m DM 


t e nd(Gyr) 


HighSoftC 


514 


3 


x 10 4 


5 


x 


10 s 


1 x 10 5 


0.148 


4.5 


MidSoftC 


200 


3 


x 10 4 


5 


X 


10 5 


1 x 10 5 


0.148 


3.9 


LowSoftC 


60 


3 


x 10 4 


5 


x 


10 5 


1 x 10 5 


0.148 


3.3 


LowSoftFloorC 


60 




300 


5 


x 


10 5 


1 x 10 5 


0.148 


3.7 


LowResC 


60 




300 


1 


X 


10 5 


1 x 10 5 


0.148 


1.6 


LowMassC 


512 


3 


x 10 4 


5 


X 


10 s 


1 x 10 5 


0.030 


7.8 



Table 3. Summary of collapse runs. l so f t is the minimum softening length, Tji oor is the temperature floor, n g and npjM are the numbers 
of gas and dark matter particles, m g /mrjM is the gas/dark matter mass ratio, and t ent j is the total simulation time. 



of B. If B contains at least half of the particles contained 
by cloud A, then B is a child of A. If B has several parents, 
then a merger has occurred. If A has several children, then 
a separation has occurred. If A is the only parent of B, and 
B is the only child of A, then B is identified as the same 
cloud as A. This categorization allows for multiple parents 
to join in a merger and it is also possible for a parent to split 
into into multiple children. During simulations we observed 
that mergers can be complex with clouds merging and sep- 
arating several times before settling into a single cloud, or 
in some cases while no longer interacting. This means our 
statistics are perhaps better thought of as recording 'inter- 
action' rates (including 'self-interaction') rather than cloud 
collision rates. 



2.3.2 Treatment of cloud energy and interactions 

To quantify the energy loss due to interactions, we compare 
the kinetic and potential energies of clouds in sequential out- 
puts across a separation or merger event. For the 'combined' 
stage of the interaction (the earlier output for a separation, 
the later output for a merger) we calculate the centre of 
mass kinetic energy, 



K, 



combined 



HE' 



i 

P 



Y 



(10) 



where P is the set of all particles 'involved' in the interaction 
(defined below), and Vi and m; are the velocity and mass 
of particle i. This is compared with the sum of the centre 



of mass kinetic energies of the clouds during the 'separated' 
stage of the interaction, 



K. 



separated 



H? m< 



> Vtli 



(11) 



where Pj is the set of all particles in cloud j and C is the set 
of all clouds involved in the interaction during the separated 
stage. 

The merged or unseparated cloud does not typically 
contain all of the particles from the clouds that formed it or 
that separated from it; there are always a number of particles 
that are expelled during the interaction, while other particles 
may accrete on to the clouds during the interaction. These 
particles carry kinetic energy, and so to ensure that we are 
measuring a real loss of energy from the system and not 
just an apparent energy loss from particles leaving the cloud 
phase, we define P by 



j 



(12) 



The final step in the energy budget calculation is to 
ensure that the energy change due to clouds moving in the 
potential well of the dark matter and baryons is accounted 
for. We take out the effects of these gravitational interac- 
tions by calculating the potential between P and all other 
particles during both dumps, and subtracting the difference 
from the kinetic energy that was calculated. 
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2.3.3 Energy analysis 

From the cloud energy budget we can obtain an estimate 
for the total time-scale for dissipation of kinetic energy from 
cloud-cloud interactions. By analogy with the star-formation 
time-scale, typically defined as tsFR = E sas /(dE*/cit), we 
define the viscous time-scale due to cloud-cloud collisions to 
be the dissipative time-scale, 



K 



(13) 



-dKc/dt' 

where K is the total rotational kinetic energy of the gas, 
and — dKc/dt is the rate at which this kinetic energy is dis- 
sipated due to collisions. This dissipation rate can be written 
as 



-dK c _ n AK col 



(14) 



where C is the interaction rate (determined by counting the 
number of interactions that occur within a time period) and 
AK co \ is the energy lost per interaction. In practice, we av- 
erage over n co i interactions so that 



— i — 



n co i J2"=i 



(15) 



where At is the time period that the n co \ interactions oc- 
curred over (and hence C = ^j^), AKi is the kinetic energy 
lost in a particular interaction i and K(ti) is the total kinetic 
energy in gas at the time of that interaction. 

It is important that we connect this method of measur- 
ing the dissipative time-scale in our models with definitions 
used elsewhere. It is commonly argued (e.g. Bell 2002) that 
the form of the viscous time-scale is 



R 



(16) 



where R is the radial coordinate and v the (effective) vis- 
cosity. 

To see how this form arises in our measurements, con- 
sider the following argument: If we neglect radial velocity, 
then the kinetic energy per unit volume of a component of 
fluid in a rotating disc is k = p(RQ.) 2 /2. We can convert 
the rate of viscous dissipation for a generic fluid (<£», the 



energy lost per unit volume per unit time) from Mihalas 



& Weibel Mihalas (19841 into cylindrical coordinates and 
again assume angular velocity dominates, simplifying it to: 



$ = pu{RQ!) 2 , 



(17) 



where the prime indicates a radial derivative. We can sub- 
stitute these values into our definition for t„ ool because 
$ = — dkc/dt, so 

_ » 2 
~~ 2u{fl') 2 ' 



p(RQ) 2 /2 
pu(RQ') 2 



(18) 



If we then take a power law for rotation Q oc R a then 

which agrees with R 2 jv within a factor of l/2a 2 . For a flat 
rotation curve, a = 1 and this factor is merely 1/2 - hence 
the dissipative time-scale is of the order of the traditional 



viscous time-scale. Note, Lin & Pringlc (1987) give a differ- 
ent prefactor - (2 — a) /(a). However, these values all agree 



within an order of magnitude, provided a is not extremely 
large or small. Although our viscous time-scales are calcu- 
lated over the whole disc to ensure sufficient numbers of in- 
teractions are measured, and the analytical R 2 jv is a local 
value at a specific radius, we should not expect this to have 
an effect beyond an order of magnitude, assuming analytical 
viscous time-scales have been calculated at a representative 
radius. 



3 RESULTS 

3.1 Milky Way Model 

3.1.1 General Evolution 

The evolution of all models excluding HighSoftMW are sim- 
ilaiQ In these models the gas disc is initially close to equi- 
librium. However, the gas rapidly cools and becomes unsta- 
ble, collapsing vertically (except in FlatMW and HighRes- 
FlatMW, which are produced from already-collapsed initial 
conditions), and forming spiral instabilities which fragment 
into large number of small (m ~ 10 6 -10 7 Mq, R ~ 100 pc) 
clouds. 

After this epoch of rapid cloud formation, the clouds 
merge and continue to accrete material. The number of 
clouds drops, as illustrated in Fig. [2] while the total mass 
within clouds continues to increase until both reach a less 
dramatic stage from around 0.8-1.0 Gyr, where the number 
of clouds decays only gradually as the mass within clouds 
gradually increases. A face-on view of the evolution of Low- 
SoftMW is shown in Fig. |3j and a snapshot of HighRes- 
FlatMW is shown in Fig. gin HighSoftMW cloud collapse 
was quenched by the high softening length, and instead the 
disc was dominated by large scale instabilities (Fig.[5|. The 
higher gas mass in HighGasMW and MedGasMW reduced 
the hydrodynamic time-step and so these simulations could 
only be run for ~ 0.45 Gyr, while the increased computa- 
tional load of the high resolution run HighResFlatMW also 
made a full simulation of 1.0 Gyr unfeasible, and so this 
simulation was evolved for ~ 0.3 Gyr. 

The gas disc separates into two phases: diffuse gas which 
retains a moderate temperature (~ 10 3 to ~ 10 4 K) though 
shock heating and a low cooling time, and dense gas whose 
temperature is tightly controlled by the Robertson-Kravstov 
dynamic temperature floor. It should be noted that while our 
models lack direct stellar feedback, the dynamic floor can 
heat the dense gas to temperatures as high as 3 x 10 4 K. This 
temperature is equivalent to a sound speed of ~ 26 km s , 
which is on the order of the velocity dispersion generated by 
various feedback mec hanisms (|Thacker fc Couchman||2000| 
Governato et al.|2007| |Ostriker fc Shetty|2011[ ). Hence while 
we expect implementing feedback would change our results, 
the difference may not be large. This is further supported 
by the findings of Shetty & Ostriker ( 2008 1 , who found 



that the properties of large clouds are not strongly sensi- 
tive to feedback. Tests were also performed with a higher 
cooling floor of 3 x 10 4 K, and no clouds were formed. This 
perhaps demonstrates that a static cooling floor is a worse 



1 Animations for some models presented here are available at 
http: // ap . smu. ca/~thacker/williams/ cloudcols .html 
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Figure 2. Number of clouds in Milky Way models. To smooth the data, each plotted point is an average of the 29 data points centred 
on it. The number peaks when many clouds are rapidly formed as the gas temperature drops below the Toomre instability threshold. It 
drops as these clouds merge. 



approximation to feedback as it inputs energy into any cool 
region of gas regardless of density, impeding any collapse 
that would have actually formed stars, in contrast with a 
dynamic temperature floor which inputs energy only into 
dense star-forming gas. 



3.1.2 Cloud formation & numerical issues 

We now draw attention to the differences between the sim- 
ulations illustrated in Fig. [2] While LowResMW produces 
clouds at the same time as LowSoftMW (top left), it pro- 
duces fewer of them as the mass spectrum is truncated. Sim- 
larly, FlatMW produces clouds at the same time as High- 
ResFlatMW, but in smaller numbers (bottom right). Hence 
there is a trend of producing more clouds with increasing 
resolution. Overall, the flat initial conditions of FlatMW 
and HighResFlatMW produced clouds earlier and in greater 
numbers than in LowSoftMW. LowViscMW appears identi- 
cal to LowSoftMW, suggesting that numerical artefacts due 
to artificial viscosity are not a significant effect (top right). 
LowFloorMW produced more clouds than LowSoftMW as 
the lower cooling floor allows the disc to become more un- 
stable to cloud formation from Toomre instabilities. We also 
found that clouds formed earlier and were more numer- 



ous with increasing gas fraction, as demonstrated by High- 
GasMW and MedGasMW (bottom left). 

We found that replacing the halo with a static potential 
did nmot have a significant effect - the mass spectra and 
number of clouds formed over 430 Myr of evolution were 
almost identical (Fig. [6]l . 

As expected, the gravitational softening parameter has 
a significant effect on cloud formation. With a softening of 60 
pc (LowSoftMW), a maximum of ~ 300 clouds were formed 
at a time of 0.02 Myr, while with a softening of 200 pc (Mid- 
SoftMW), half as many were formed (~ 150), and the peak 
number was achieved later (0.04Myr). It should be noted 
though, that both models have a similar fraction of mass 
in clouds (~ 80%). Increasing the softening yet further to 
500 pc (HighSoftMW), leads to almost no clouds forming 
other than a few clouds in the centre of the galaxy after 
about a Gyr of evolution (not shown in Fig. These re- 
sults match what would be expected on theoretical grounds. 
Increasing the softening length delays cloud formation and 
produces fewer, more massive clouds, unless the softening 
length is increased above a certain threshold, beyond which 
cloud formation is prevented. 

It seems most likely that this threshold softening length 
is related to the wavelength of the unstable mode that causes 
cloud formation. We can calculate this using the two-fluid 
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Figure 3. Evolution of LowSoftMW. A featureless disc (top-left) rapidly collapses into a larger number of clouds (top-right) after 
around 200 Myr of evolution. These clouds interact with each other and accrete material from 400 Myr (bottom-left) until the simulation 
ends after 1.1 Gyr (bottom-right). 



(gas/star) Q gs stability parameter from Jog & Solomon 
(19841), iRafikovl d2001b, and |Li et al.l d2005b . The individ- 



ual Q parameters for stars and gas are denned as 



Qs = 



7rGE s 



7tGE e 



(20) 



where E s and E g are the gas and stellar surface densities, a s 
the stellar radial velocity dispersion, c g the gas sound speed, 
and k the epicyclic parameter. Note that Q s differs from 



Toomre ( 1964 1 's definition of Q for a collisionless system by 



a factor of 3.36/7T. If we define 
q = 2na s /(K\i),f = Cg/a 3 , 



(21) 



where Xi is the wavelength of a particular mode of instabil- 
ity, and treat the stars as a fluid with sound speed equal to 
a s as in Rafikov ( 2001| ) (who follows [Jog fc Solomon|1984[ |, 
we can define a combined Q gs by 



1 

Qgs 



fq 



1 + q 2 



1 + g 2 / 2 



(22) 



with a stability condition of Q gs < 1. 

We calculate Q g3 by using azimuthal means of Q, E, 
k, c g and a a, and setting Xi to X m in, the wavelength that 
minimizes Q ga - It is worth cautioning that these parameters 
are derived from linear perturbation theory and may not 
adequately describe the system once clouds have formed. 
Nevertheless, Xmin does not rapidly vary from t = 1 Myr to 
t — 200 Myr for LowSoftC as shown in Fig. [7] A m ; n is fairly 
small (< lfcpc) until a radius of 10 kpc at which point it 
triples in size. This jump is due to the small wavelength gas 
instabilities starting to dominate over the large wavelength 
stellar instabilities. A comparison with the face-on density 
plots (e.g. Fig. [3| shows that clouds form within 10 kpc. In 
this region X m i n is of the order of 100s of pc. The 'threshold' 
resolution for cloud formation (assuming 4 to 5 softening 



© 2011 RAS, MNRAS 000,[T|fT8l 



Effective viscosity from cloud-cloud collisions 

Static Potential 




Mass (Solar masses) 



in three-dimensional global SPH simulations 9 



Static Potential 




0.1 0.2 0.3 0.4 0.5 

Time (Gyr) 
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Figure 4. HighResFlatMW after 300 Myr of evolution 



lengths are required) in our models lies somewhere between 
200 pc and 500 pc, and is consistent with this range. This 
quantifies an often quoted caveat for galaxy models - if the 
gravitational softening length is larger than the wavelength 
of the most unstable modes, then fragmentation is artificially 
frustrated. 

The size of the unstable perturbations can be used to 
crudely estimate the masses of clouds. Assuming that the 
disc fragments into clumps of mass ~ 7rEA^j n , then for 
the LowSoftMW simulation (for example) the typical cloud 
masses should be the order of several 10 6 Mq, which is ad- 
mittedly significantly larger than average molecular cloud 
masses and actually much closer to giant molecular cloud 
complex masses. Nonetheless, this value is broadly consis- 
tent with our spectrum of cloud masses (e.g. Fig. [8}. How- 
ever, we caution against over interpretation as the mass 
spectrum convolves together an initial spectrum and its sub- 
sequent evolution. If this simple approach to calculating ini- 



Figure 5. HighSoftMW after ~ 1.5 Gyr of evolution. Because of 
the large softening length, the disc does not undergo local frag- 
mentation into clouds, and is instead dominated by bar and spiral 
instabilities. 



tial cloud masses were accurate we would not expect a higher 
resolution model to produce smaller clouds from this mode 
of instability, although non-azimuthally symmetric modes 
which may produce smaller scale instabilities have been ex- 
cluded from this analysis. Smaller clouds could also be pro- 
duced in a higher resolution Milky Way model by changing 
the initial conditions, or if these giant clouds undergo further 
fragmentation. 



3.1.3 Cloud Mass Functions 

The mass functions of our clouds (Fig. 



of Tasker fc Tan (20091 and Agertz et aid 1 2009 1 in that our 



differ from those 



clouds are more massive. However, neither of these studies 
has equivalent physics. Tasker & Tan differ in that they do 
not include dynamic stellar disc while Agertz et al include 
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feedback. Resolution could potentially also be an issue: al- 
though our mass function does not greatly vary between our 
low and moderate resolution models in our fiducial simula- 
tions, our high resolution flat model produced lower mass 



clouds than the moderate resolution flat model (Fig. 10 1. 

The high-mass region of our cumulative mass spectrum 
plot follows a power law (i.e. N(m) oc m a or N(m > 
M) oc M a+1 ). A least-squares fitting gives a ~ —1.5. This 
is slightly shallower than the ~ —1.8 in the simulations of 
Das fc Jog| {T996| >, and |Dobbs fc Bonnell| pOOS) but close 
to the values of —1.5 to —1.6 from observations (Sanders 
et al.||1985| |Solomon et aT][Tg87| [Solomon fc Rivolol|1989| 



Williams & McKce 1997 



Roman- Duval et al. 



2010), and 



from the simulated mass spectra at around 10 Mq at 300 



Myr in Tasker fc Tan ( 2009 1 and at 1 Gyr in Agertz et al 



(2009) 



3.1.4 Viscous time-scales 

The viscous time-scale is calculated using the method de- 
scribed in section |2.3| and is plotted in Fig. |11| Each point 
is calculated from 600 collisions. There is a general trend 
toward lower time-scales as the simulation evolves, and the 
final time-scales are generally below 10 Gyr, with many ap- 
proaching 1 Gyr. This decreasing trend coincides with a 
trend of the number of clouds lowering and the mass of in- 
dividual clouds increasing. The time-scales are less than a 
Hubble Time, and so should have some significant effect on 
the evolution of a galaxy, contrary to the predictions of B02. 

This energy loss is seen in the mean specific kinetic en- 
ergy of the gas in LowSoftMW, dropping from 1.9 x 10 14 
erg/g at t = 170 Myr to 1.1 x 10 14 erg/g at t = 1010 Myr. 
As expected, this loss is primarily in the clouds - the diffuse 
gas only drops from 1.9 x 10 14 erg/g to 1.7 x 10 14 erg/g in the 
same time period. We can make an additional crude estimate 
of the total viscous time-scale, t v = {t% — ti)ki/(ki — fe), 
where ti and fc, are the time and specific kinetic energy at 
these two outputs. This gives a time-scale of t v ~ 2 Gyr, 
around half of the value of t v = 4.5 Gyr from our method in 
section [273] This is either because additional energy is being 
dissipated through internal processes in clouds, or because 
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Table 4. Mean viscous time-scales and simulation lengths for all 
runs for the time from the first to the last recorded interaction. 
These time-scales are the mean time-scales during the time pe- 
riod from the first to the last recorded interaction, time-scales are 
not given for LowMassC and HighSoftMW. There were no inter- 
actions in LowMassC, as it did not form clouds. Interactions were 
detected in HighSoftMW, but only in clumps within the central 
bar, which do not contribute to disc viscosity. The viscous time- 
scale for the first 300 Myr of LowSoftMW and FlatMW are also 
given for more direct comparison with HighResFlatMW. 



interactions are being missed by our interaction finding pro- 
cedure. The energy loss of the gas over this time period in 
LowViscMW is the same to a precision of less than 0.1 x 10 14 
erg/g, suggesting that this additional energy loss is 'physi- 
cal' and not dominated by artificial viscosity. Nevertheless, 
we can not evaluate how much of this additional energy loss 
is directly due to cloud-cloud collisions, and so it is more in- 
formative to use the procedure in |2.3| to calculate the viscous 
time-scale. 

The mean viscous time-scales from all interactions over 
each entire simulation for both the Milky Way and collapse 
models are tabulated in Table|4] Despite the variation of pa- 
rameters, many of the time-scales are within a narrow range, 
from 3-5 Gyr. Modifying the artificial viscosity (LowVis- 
cMW) did not appear to significantly change the viscous 
time-scale. The softening length in HighSoftMW (600 pc) 
was large enough to completely quench cloud formation, ex- 
cept for a few clumps that formed within the central bar in- 
stability. We do not include a viscous time-scale here as the 
mechanisms for formation and interaction are different to 
those of molecular clouds in nearly circular orbits. Feedback 
processes from star formation and AGN would also be more 
important here than in the other models. However, lowering 
the softening length from 100 pc to 60 pc (MedSoftMW to 
LowSoftMW), while increasing the number of clouds pro- 
duced, did not significantly alter the viscous time-scale. 

HighGasMW has a significantly shorter viscous time- 
scale at 0.6 Gyr, and indeed there appears to be a trend of 
decreasing viscous time-scale with increasing gas fraction. 
This is clearer if we compare the models over the same time 
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Figure 8. Mass spectra for clouds in Milky Way runs at 800 Myr. Left: Cumulative mass spectra (for comparison with |Agertz et ah] 
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Figure 9. Cumulative cloud mass spectra across runs with vary- 
ing gas fraction. 



Figure 10. Cumulative cloud mass spectra from flat initial con- 
ditions, including our highest resolution model. 



period. The viscous time-scale over the first 430 Myr is 7.1 
Gyr for LowSoftMW, 1.5 Gyr for MedGasMW, and 0.6 Gyr 
for HighGasMW. Increasing the gas fraction increases the 
mass of the cloud population (Fig. |9j, which increases the 
frequency and dissipative efficiency of collisions. 

HighResFlatMW is our highest resolution simulation, 
but has different initial conditions to LowSoftMW due to the 
more stringent stability requirements at high resolution (de- 
tailed in section 12. 2.1}. The flat discs of FlatMW and High- 



ResFlatMW caused cloud formation to occur earlier than in 
LowSoftMW. A resolution dependence is also evident: The 
2.5x increase in mass resolution from FlatMW to HighRes- 
FlatMW caused a 1.4x increase in the viscous time-scale, 
and the 5x increase in mass resolution from LowResMW to 
LowSoftMW caused a 1.8x increase in the viscous time-scale. 



3.2 Monolithic Collapse Model 

In all models the gas collapse proceeds as soon as cooling is 
turned on, thus breaking the hydrostatic equilibrium. The 
High-S profile slowed the collapse sufficiently for the infalling 
gas to fragment into clouds at a large radius, although these 
clouds are are too diffuse to be found by the cloud identifi- 
cation algorithm. As the simulation progresses, these clouds 
start to merge (from t ~ 3 Gyr in all runs except for Low- 
MassC), and reach the effective threshold density of our 
cloud-finder. The number of clouds quickly reaches a maxi- 



mum (see Fig. 15 1. These clouds combine to form a disc. The 
number and size of clouds these discs fragment into varies 
greatly between our models. 

In HighSoftC, MidSoftC, LowSoftC and LowResC, the 
disc is extremely unstable, collapsing int o ~ 7 massive (sev- 
eral times 10 9 Mq in mass) clumps (Fig. 121. These are not 



small-scale GMC-style clumps as found in the Milky Way 
simulations, and perhaps this level of collapse is more anal- 



© 2011 RAS, MNRAS 000,[T}{T8] 



12 D. J. Williamson and R. J. Thacker 



10' 



8 3 io 10 

5 5> 



10 M 



□ 



0.2 



LowSoftMW — h- 

MedSoftMW —x- 

LowFloorMW x- 

LowViscMW e 

LowResMW — -»- 



0.4 0.6 
Time (Gyr) 



-a 



0.8 



10 



10 



8 8 io 10 



10 M 



10° 



Q 

H 




MedGasMW - 
HighGasMW 
FlatMW 
HighResFlatMW 


— I 

— -X— - 

---X 

H~ 


\ % 










^ a 










X _ 






""-x''' 







0.1 0.2 0.3 0.4 0.5 0.6 
Time (Gyr) 



Figure 11. Viscous time-scales for disc models that ran for > 800 Myr (left) and ^ 800 Myr (right). At early times, some models give 
negative time-scales, but as these values are large, they are not as dynamically important and are not plotted. 



ogous to the gas-rich clump-cluster galaxies found at high 
redshift (e.g. Elmegreen fc Elmegreen|2005 1. In the simula- 



tions of Bournaud et al. (20071 andlDekel et al. (20091, the 



large clumps in clump-cluster galaxies coalesce into a central 
bulge, forming a more stable disc. These simulations differ 
to ours particularly in that they include star-formation and 



feedback. With infalling material, Dekel et al. (2009) finds 



the clumpy phase can last for several Gyr. 

The heavy clustering in this discs dictated that they 
could only be evolved for < 1 Gyr after formation (which 
takes ~ 3 Gyr) due to problems with the SPH solver. The 
high densities cause a large increase in the number of parti- 
cles with smoothing lengths at the minimum allowed which 
contributes to an 0(n 2 ) slowdown. 



The simulations of Kaufmann et al. (20091, while in- 



cluding star-formation (but not explicit feedback) , also pro- 
duce a disc with large-scale gravitational instabilites. Both 
our and Kaufmann's models have a temperature floor of 
3 x 10 4 K, as a very crude form of feedback. Including star- 
formation and more self-consistent feedback method could 
produce a stable disc ( Stinson et al.|2006 Christensen et al 
2010). 



In LowSoftFloorC, the low temperature floor allows the 
halo clouds to condense into dense (n ~ 10 4 — 10 5 cm~ 3 ) 
clumps (Fig. 13 1. Their low cross-section means that their 
coalescence has properties of a collisionless collapse. So in 
addition to an unstable disc, there exists a swarm of clumps 
with a half-mass height of 7.8 kpc. Their ellipsoidal distribu- 
tion and high densities are reminiscent of globular clusters, 
but the inclusion of feedback would definitely increase the 
cloud cross-sections and produce a more dissipated and flat- 
tened disc. 

LowMassC is the only run that produces a disc that 
does not collapse into large clumps (Fig. 
took considerably longer to form (~ 4.5 Gyr) 



141, although it 
and the disc 

is still dominated by spiral instabilities. Discs are unstable 
to bar formation when the disc mass fraction is greater than 
the spin parameter (rrid > Ag) ( jEfstathiou et al.|1982 Foyle 
et al.|[2"o"08| ), so a lower mass disc is more stable. If the bar 
is too strong, it may fragment into large clumps. This in- 
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Figure 14. LowMassC at t = 6.0 Gyr. The disc undergoes spi- 
ral instabilities but does not fragment into clumps as the other 
collapse models do. 



stability may well drive the infalling clouds into a few large 
clumps in the higher mass models. 

As seen in Table |4) the viscous time-scales for the col- 
lapse runs trend toward lower values than the Milky Way 
simulations - around 1 — 2 Gyr. Though the number of in- 
teractions is not as large as in the Milky Way models, they 
occur over a short period (e.g. all 566 interactions in Low- 
SoftC are within ~ 500 Myr). The number of clouds is small, 
so each cloud undergoes many collisions, producing a short 
viscous time-scale. 



3.3 Comparison with Analytical Model 

B02 argued that while cloud collisions are not uncommon 
(occurring £ 1 time per orbit), the low efficiency of cloud 
collisions produces a long viscous time-scale. This efficiency 
is measured with a parameter rj, equal to the fraction of a 
cloud's energy that is lost in a collision (not entirely dis- 
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Figure 12. Impact of varying the softening length and resolution in collapse runs at t=3.5 Gyr. Top left is HighSoftC (514 pc, 3 X 10 4 
K), top right is MidSoftC (200 pc, 3 X 10 4 K), bottom left is LowSoftC (60 pc, 3 X 10 4 K) and bottom right is LowResC (60 pc, 300 
K). Although HighSoftC, MidSoftC and LowSoftC produce different numbers of clouds initially (more clouds for a shorter softening 
length), after ~ 500 Myr of collisions all three models have ~ 7 large clumps. Despite the low temperature floor, the limited resolution 
of LowResC produces an unstable disc, instead of a swam of dense clumps as in LowSoftFloorC. 



similar from a coefficient of restitution). When two clouds 
merge completely, the fraction of kinetic energy lost is well 
approximated by r\ — (v TC i/v TOt ) 2 , where t; rc i is the relative 
velocity of the clouds, and v TO t is their rotational velocity 
which is roughly constant for a galaxy. This is consistent 
with our numerical results. The analytical model of B02, 
finds that r\ < 10~ 2 for a Milky- Way-like model, concluding 
that cloud-cloud collisions are not an efficient sink of energy, 
with ty ~ 1000-2000 Gyr. 

The complex interactions that occur between clouds 
in our simulation mean that it is not straightforward to 
determine values for r\. Several of our merger and sepa- 
ration events can take place within what is really a sin- 
gle extended interaction, which lowers the average time be- 
tween interactions significantly. Indeed, we find the inter- 
action rates are on the order of one separation or merger 



event per cloud every 50 — 60 Myr for LowSoftMW, Med- 
SoftMW, LowFloorMW, LowViscMW, HighSoftC and Low- 
MassC. The greatest interaction time-scale was in LowSoft- 
FloorC (335 Myr), and the smallest was in LowSoftC (14 
Myr). 

It is difficult to track the number of interactions over a 
full merger process, as additional clouds often interact with 
the merging clouds. We carefully examined examined a span 
of time around each of a sample of 10 recorded interactions in 
LowSoftMW on an iteration-by-iteration basis to determine 
the number of recorded interactions per 'real' interaction. 
These interactions were selected so that they were evenly 
distributed across the simulation (~ 2 every 5000 iterations). 
We initially examined a period of ±800 iterations around the 
interaction, and if no 'real' interaction was observed during 
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Figure 13. Face-on and side-on density plots of LowSoftFloorC at t=3.7 Gyr. The swarm of clumps has a half-mass height of 7.8 kpc. 
The disc is very chaotic: 10 kpc, the azimuthally averaged tangential velocities and velocity dispersions are 180 km s- 1 and 105 km s^ 1 . 
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Figure 15. Number of clouds in collapse models (excepting Low- 
SoftFloorC). To smooth the data, each plotted point is an average 
of the 29 data points centred on it. Being very unstable, these sys- 
tems formed a few large clumps rather than many small clumps. 



this time, this was extended to ±2000 iterations. Several 
different behaviours were observed: 

• In two cases, no real interactions were observed; outer 
parts of a cloud were attaching and detaching to the main 
cloud, and dissolving and condensing across the cloud den- 
sity threshold, causing a number of recorded interactions 
which did not correspond to any clear long-term merger, 
scattering or separation event. 

• Three events were 'messy' interactions with 6, 7 and 16 
recorded interactions per real event; the event consisting of 
16 recorded split and merge events was a scattering event 
where the clouds passing by each other several times before 
separating for a final time. 

• Four more events were more 'tidy' interactions, with 1, 
2, 3, and 4 interactions per real event. 



• The last event was a series of mergers in rapid succes- 
sion - 3 recorded and 3 'real' mergers. 

Overall, there was a mean of 4.9 recorded interactions 
per examined period, with a standard deviation of 4.3. 
A total of 11 'real' interactions were observed, giving 4.5 
recorded interactions per real interaction. This increases our 
interaction time-scale to one event per ~ 250 Myr for the 
LowSoftMW-like models. This is approximately once per or- 
bit at a solar radius. The analytic estimate in B02 of the 
cloud-cloud collision rate is ~ 100 Myr, which is of similar 
order. 

We can estimate an rj for the interactions in our models 

by 



V = -(AK + A<f>)/(K C ), 



(23) 



where AK and A<jf> are the change in kinetic and poten- 
tial energy of a cloud, and K c is the total kinetic energy of 
both clouds before collision. rj can be negative, as energy is 
converted from internal motions into orbital kinetic energy 
during separations. The clouds all have similar velocity be- 
cause of the flat rotation curve, so the total energy lost is 
primarily dependent on rj and the cloud masses. We find for 
most interactions \rj\ is on the order of ~ 0.002 (Fig. 16 1. If 
we separate our rj values into two sets, rj- for rj < and rj+ 
for rj > 0, we find that the median value of \ij-\ is greater 
than the median value of 1 77+ 1 , even though the viscous time- 
scale is positive. This is because although rj, the relative en- 
ergy change is larger for interactions which increase orbital 
energy (rj- < 0) than those which decrease orbital energy 
(77+ > 0), the absolute change in energy is larger for inter- 
actions which decrease orbital energy than increase it - i.e. 
these interactions tend to occur between clouds with greater 
mass. Although it is not apparent on these plots, there are 
several collisions for which rj is very large, with rj > 0.1. 
These interactions occurred within the 1 kpc of galaxy cen- 
tre, and only after ~ 400 Myr. These are clouds that have 
been strongly scattered by interactions and fallen down the 
potential well, colliding with speeds of > 100 km s _1 . 
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Figure 16. Distributions of the fraction of energy lost in a collision r), in bins of 0.35 dex. Left: LowSoftMW, Right: LowSoftC. For 
each simulation, the distribution of all r\— < and ?y+ > are plotted separately. In both cases, the median value of \r\— | is greater than 
the median value of |r;+|, even though both models show a positive viscous time-scale. 



Our interactions are no more efficient at removing en- 
ergy than in B02, and are no more common, yet the B02 
model predicts t v ~ 1000-2000 Gyr, while our simulations 
have t v < 10 Gyr. Our simulated discs are more energetic 
than standard Milky Way models: the velocity dispersion 
in LowSoftMW is ~ 20 km s _1 at 7.5 kpc, more than triple 
the standard Milky Way value used in B02 (6 km s _1 ). How- 
ever, this is not the cause of the large difference between the 
model of B02 and our own. Here we derive our own model 
for 77, and contrast this with the model in B02 to find the 
source of this disparity. 

We can split the velocity components of v re i into tan- 
gential and radial components to give 



Viel 
Viot 



R 2 ((f>i - (P2) 2 + (Rl - R2) 



(24) 



If we make the epicyclic approximation) Binney & 
Tremainc 2008), that the deviation from a circular orbit is 
small compared to the radius of the orbit (7? = Rg+x, where 
R s is the 'guiding centre' of the orbit, and x <C R is the ra- 
dial excursion), then R — x = Xncos(Kt + a), (where X is 
the maximum radial excursion of a cloud, kappa the epicyclic 
frequency, and a is a phase parameter) and <j> — R g v TOt /R 2 
from conservation of momentum in a flat rotation curve. 
Hence R 2 (if>\ — 4>2) 2 = (v 2 ot /R 2 )(R g ,i — R s ,2) 2 - the tangen- 
tial component of the difference in velocity depends only on 
the radial distance between the clouds' guiding radii. 

The radial component is more difficult to calculate, as 
it depends on the phase of the interaction. We can estimate 
the maximum rj by assuming the clouds are perfectly out of 
phase - that is, 



(ill 



R2) 2 ~ (XlKl + X2K2) 2 



2vi 



2vt 



R 2 



Xi X2 
R s ,i R g ,2 



-(X1+X2), (25) 



as R ~ -Rg. For clouds to collide precisely out of phase, they 
must have the same guiding radius, and so _R gi i — i? gi 2 = 0. 
Hence, if Xi ~ X 2 ~ X, then ?y max , r = 8X 2 /R 2 . If the clouds 
are at their maximum deviation when they collide, then their 
radial velocities are zero, but their relative cj> velocities are 



maximised, that is, 4>\ — 4> 2 = 2X, and so Tj m ax,ij> = 4X 2 /7i 2 . 
These coefficients give the maximum rj, but we should nev- 
ertheless expect rj ~ X 2 /R 2 , i.e. rj depends on the radial 
excursion of clouds. 

This can also be expressed in terms of a velocity disper- 
sion. We can calculate the velocity dispersion by 



<(v- 



,) 2 ) = m 2 ) + (R 2 (<b-n g ) 2 ) 



(26) 



Assuming a flat rotation curve and that X and ft are 
more or less constant within the region of interest, the radial 
component is x = Xncos(Kt + a) , hence 



<(±) 2 > = (1/2)XV =X 2 v 2 ot /R 2 , 



(27) 



and the tangential component is R(<fr — fl g ) — 
—2XQ, g sin(fd; + ct), hence 



{R 2 {j> - n g ) 2 ) = 2X 2 fl g = 2X 2 v 2 ot /R 2 . 

This gives 
vl = 3v 2 ot (X 2 /R 2 ), 
and so 



(28) 



(29) 



(30) 



From these expressions for rj we can determine the dis- 
sipative time-scale from t v = t c /r). 

We next summarise the model of B02. In the limit of 
rapid collisions, the kinematic viscosity due to cloud-cloud 
collisions can be modelled as a Reynolds stress and expressed 
as v ~ \dv a ( |Faber|1995[ ), where v s is the velocity dispersion, 
and Xd is the mean free path. The mean free path is Ad = 
v a t c , where t c is the typical time between collisions. Similarly 
to our result, B02 states r\ ~ AR 2 /R 2 , where A_R is the 
radial distance between collisions. For the case of very rapid 
collisions, AT? ~ Ad, so rj ~ \ 2 d /R 2 . This gives a viscous 
time-scale that should be equal to the dissipative time-scale, 



tu 



IF 

V 



R 2 

XdV s 



■ \ 2 



(31) 



Hence if we follow the description given in B02, the 
results should be equivalent to ours. Continuing to follow 
B02, we can set t c = M c i oud h/(T, g v a irr 2 loud ), so 
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R 2 



R 2 7rr cloud^'g 



(32) 



Fit 



We can evaluate this using the Milky Way parameters 
of B02, that r = 7.5 kpc, v TO t = 220 km s" 1 , v s = 6 km 
s" 1 , E g = 5OM pc" 2 , A'/ c i oud = 1O 5 M , h = 100 pc, and 
r cloud = 10 pc to result in t v = 14 Gyr. However, B02 states 
t v ~ 2000 Gyr. This disagrees by a factor of 1 /r). It appears 
that B02 includes an additional factor of r\ ~ 0.008 in the 
denominator - i.e. t v ,Beil ~ -R 2 /(V u )- This 77 is not necessary, 
as it is already included in the radial excursion or velocity 
dispersion, and as is clear from equation |31[ the expression 
t v ~ R 2 /vq is not equivalent to the dissipative time-scale. 

In B02 's rare collision case, v ~ v 3 AR(t K /t c ), where 
t« = 27t/k is the epicyclic time-scale. For a flat rotation 
curve k = V2Q ~ v /i?. The excursion Ai? is on the order 
of the radial excursion of the epicyclic motion of the clouds. 
B02 state AR ~ v a /n ~ v s R/v Iot , which is consistent with 
our result in equation |29| Putting this together gives 



t c /(2w)v 



2 

rot 



(33) 



i.e. r\ ~ 27ru 2 /(i) 2 ot ) ~ 0.023. B02 uses a low surface 
brightness galaxy in this case, with E g = 1OM pc -2 and 
«rot = 100 km s _1 , which results in t v ~ 23 Gyr. Again, 
the value in B02 is much larger, t„ ~ 1000 Gyr, which again 
is higher than our calculated value by a factor of approxi- 
mately I/77. 

These model are intended to apply in the limits of very 
frequent or very infrequent collisions where t<Sl < 1 or 
tcQ. 1. In our simulations, we found that clouds collide 
about once per orbit, i.e. Qt c ~ 1. However, we can contrast 
these results with those of |Goldreich fc Tremaine| ( |1978[ ) , 
who solve the Boltzmann equation for a system of inelasti- 
cally colliding particles in a disc, and find for arbitrary Qt c 
that the viscosity is of order 



v 3 \d 



1 + (fii c ) 2 



(34) 



- 20 km s _ \ 

35 pc, and Af c i oud ~ lO 7 A'/ ) for LowSoftMW at 



after we make the substitution that Xd ~ v s t c . For Qt c = 1, 
v = l/2(Ad« s ). The frequent collision case of B02, v ~ v s Xd, 
is accurate to this within an order of magnitude if we exclude 
the erroneous factor of I/77. 

Substituting our typical cloud and disc parameters at 
7.5 kpc (h ~ 25 pc, E g ~ 100M Q /pc~ 2 

^"cloud 

t = 1 Gyr into this model gives a viscous time-scale of 1.1 
Gyr. This values somewhat underestimates our numerical 
results for the Milky Way models in Table [4] for most of 
which t v 1? 4.0 Gyr. The unstable disc of LowSoftC, forming 
from a collapse without stars, has very different properties 
at R = 7.5 kpc, with h ~ 250 pc, E g ~ 5OOOM /pc -2 , 
v s ~ 100 km s _1 , r c i oud ~ 100 pc, and Af c i oud ~ 1O 9 M . 
This gives t v — 0.35 Gyr, which agrees with our simula- 
tion result (0.8 Gyr) within a factor of ~ 2. The analytical 
expression for i„ was evaluated from order-of-magnitude ar- 
guments and assumptions that may not be entirely valid 
in our simulations - particularly in models with very few 
clouds, such as LowSoftC. Numerical factors also vary our 
simulation results by a factor of ~ 4. Given these issues, it 
is not surprising that the agreement is not exact. 

Interestingly, despite the different disc properties, Low- 
SoftC and LowSoftMW have similar viscous time-scales in 
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Figure 17. Correlation between peak number of clouds 
(^Vcloud max) an d viscous time-scales (t v ) for all models whose 
time-scale is given in Table [4] The plotted fit is t„ = 
(0.67 Gyr)(7V cloud , max ) - 39 . 



both our numerical simulations and in this analysis. This 
is to be expected from equation |32| We should expect the 
typical cloud mass to increase with the gas density and typ- 
ical cloud radius, and so M c i oud /(7rr 2 loud E g ) should vary 
only weakly. Hence the viscous time-scale will primarily 
depend primarily on h and v s . This suggests that time- 
scales will not vary greatly for models beyond those sim- 
ulated here - perhaps even of higher resolution. To quan- 
tify this, we note that there appears to be a correlation be- 
tween the maximum number of clouds formed (iV c i u d . max ) 
and the viscous time-scale (Fig. 
a power-law t v oc (iY c ioud,ma X ) m 



17) 



Performing a fit to 
we find a power index 
of m — 0.39 ± 0.19. This predicts a viscous time-scale of 
t v ~ 23 Gyr for iVcioud.max = 10 4 , and t v ~ 60 Gyr for 
Afckmd.max — 10°! although we caution that this is a purely 
empirical fitting and is not likely to be very accurate. 



4 CONCLUSIONS 

Previous estimates of the viscous time-scale suggest that 
the viscous time-scale for cloud-cloud collisions in a Milky- 
Way-like galaxy is large, with t u > 1000 Gyr. To test the 
hypothesis that the viscous time-scale is long, we performed 
simulations using the AP 3 M SPH code HYDRA with cooling 
down to 10 K and a dynamic temperature floor. The simu- 
lations fell into two sets of models: initially stable gaseous 
discs within a dark-matter haloes and stellar discs, and a 
gaseous spheres collapsing inside dark-matter haloes. These 
two models were chosen to bracket a wide range of stabil- 
ity. The viscous time-scale was measured by tracking clouds 
with a friends-of-friends algorithm, and determining the en- 
ergy loss when clouds collided. 

Although our cloud masses are larger than those found 
in other simulations, potentially due to insufficient resolu- 
tion, a simple analysis suggests that we are resolving the 
wavelength of the most unstable mode. However, further in- 
stabilities (in particular, non-axisymmetric instabilities that 
we exclude from our analysis) may appear at higher resolu- 
tions, and while the inclusion of energy input from stellar 
feedback may not greatly alter the properties of clouds, it 
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may contribute to cloud evaporation and affect their colli- 
sional behaviour by increasing their cross-section through 
heating. 

Identifying clouds and interactions between clouds is 
still a difficult task, as clouds have complex structures and 
dynamics. The friends-of-friends algorithm often identifies 
clouds as merging and separating several times over a pe- 
riod that upon visual inspection appears to be a single in- 
teraction. Through a detailed examination of 10 interac- 
tion events, we determined that each 'real' interaction corre- 
sponds to ~ 4.5 interactions found by our algorithms. This 
also complicated our estimates for 77 = AK c \oud/K c \o\id, the 
efficiency of energy loss per cloud interaction. We found that 
despite our low viscous time-scales, n was not large, with 
77 ~ 0.002 per recorded interaction. 

Most models from both sets of initial conditions col- 
lapsed into discs dominated by clumps of gas. The Milky 
Way models produced a more stable disc with a large num- 
ber of small clouds, while the collapse models produced a 
highly unstable disc consisting of a small number of mas- 
sive clumps. Despite this large disparity, the viscous time- 
scales were similar, with t v = 4.5 Gyr for LowSoftMW, and 
t v = 0.8 Gyr for LowSoftC. These values are much smaller 
than estimates using the formulation of B02, which overes- 
timate the viscous time-scale by appearing to erroneously 
include inefficiency of cloud collisions twice. Removing this 
factor gives analytic estimates of t v — 1.1 Gyr for Low- 
SoftMW and t v = 17 Gyr for LowSoftC. These values do 
not exactly coincide with our measured values as they are 
based on simple arguments that are particularly inaccurate 
for LowSoftC. However, they all agree with the general state- 
ment that viscosity due to cloud-cloud collisions is not neg- 
ligible. 

The scatter of t v across our models was quite small 
(0.6-16.0 Gyr), despite the range of cloud properties. Hence 
our viscous time-scales are applicable for a wider range of 
galaxies than those modelled here, although viscous time- 
scales will likely increase somewhat as resolution improves. 
For a simulation capable of resolving 10 clouds, we predict 
a viscous time-scale of around 60 Gyr, admittedly making 
the effect comparatively weak within a Hubble time, but 
nonetheless over an order of magnitude faster than previous 
estimates. 

These results suggest that viscosity due to cloud-cloud 
collisions, while not dominant, does not have a completely 
negligible effect on the evolution of a galaxy. Although our 
models may underestimate the viscous time-scales due to 
resolution effects, it still appears that cloud-cloud viscosity is 
more significant than previously estimated. While numerical 
models of galaxies may be able to model this directly (as we 
do in this work), it may be necessary to include a cloud- 
cloud viscous term in analytical and semi-analytical models 
of disc evolution. 
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